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1 Overview 


This report summarizes a project sponsored by the NASA Langley Research 
Center through the Langley Grad-Aero Program under Grant No. NAG-1-1509. 
The project was approved for the period May 6, 1993 - November 4, 1997. However, 
due to reductions in NASA’s budget and limited funds for the Grad- Aero Program, 
funding was halted in November 1994. A one-year no-cost extension was requested 
so that the student could attend the AIAA 1995 Aerospace Sciences Meeting. Dur- 
ing that period supplemental funding was made available to support the student 
over the summer of 1995. The extended termination date for the NASA grant was 
November 5, 1995. The project has been continued by supporting the student as a 
Teaching Assistant within the School of Aeronautics and Astronautics. Although 
funding for this project was terminated, NASA Langley has continued to provide 
access to their computer facilities. 

The objective of this project was to evaluate and develop subgrid-scale (SGS) 
turbulence models for large eddy simulations (LES) of compressible flows. During 
the first phase of the project results from LES using the dynamic SGS model 
were compared to those of direct numerical simulations (DNS) of compressible 
homogeneous turbulence. The findings were published in Ref. [1-3]. Ref. [3] is 
included in this report as Attachment A. 

It is becoming apparent within the LES community that numerical errors can 
have a significant impact on large eddy simulations. Ref. [4] reports results of a 
study on the effect of the formulation of the governing equations on aliasing errors. 
A copy of the galley proof for this paper is included as Attachment B. 

The second phase of the project involved implementing the dynamic SGS model 
in a NASA code for simulating supersonic flow over a flat-plate. The model has 
been successfully coded and a series of simulations has been completed. One of 
the major findings of the work is that numerical errors associated with the finite 
differencing scheme used in the code can overwhelm the SGS model and adversely 
affect the LES results. A brief write-up of the results [5] has been submitted to the 
AIAA 1997 Aerospace Sciences Meeting and is included as Attachment C. 

One of the goals of the Langley Grad- Aero Program is to train young researchers 
in Aerospace Engineering and related fields. The student supported by this project, 
Evangelos T. Spyropoulos, has performed very well. He has grown as a researcher in 
the very difficult field of turbulence. In 1993 he participated in the ICASE/NASA 
Langley summer workshop on Transition. Turbulence and Combustion. He again 
visited NASA Langley during the summers of 1994 and 1995. It has been very 
beneficial for him to interact with NASA personnel and with NASA visitors from 
all over the world. Evangelos is currently finishing his Ph.D. dissertation and should 
graduate this August. He is looking forward to a career in the aerospace industry 
where he can put his understanding to good use. 
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Evaluation of the Dynamic Model for Simulations of 
Compressible Decaying Isotropic Turbulence 


Evangelos T. Spyropoulos* and Gregory A. Blaisdelf 
Purdue University. West Lafayette . Indiana 47907 


Several issues involving the use of the dynamic subgrid-scale model in large-eddy simulations of compressible 
turbulent flows are investigated. The model is employed in simulations of compressible decaying isotropic tur- 
bulence, and its performance is compared against results from direct numerical simulations and experiments. 
Results from a parametric study suggest the model captures compressibility effects well. Use of the dynamic model 
in simulations of inhomogeneous flows requires filtering of the flowfleid In physical space rather than Fourier wave 
space. The use of spatial filters is examined by conducting simulations of isotropic turbulence. Several implicit 
filters are found to perform extremely well and similar to the sharp cutoff Alter. One explicit filter performed well, 
but all others provided excessive dissipation at higher modes. Two formulations of the dynamic model, proposed 
by Moin et ai. and Lilly, perform well, with Lilly’s being more accurate. Results suggest also a great Insensitivity of 
the model on the filter width ratio. A modification of the convective terms in the momentum and energy equations 
is found to reduce the effects of aliasing errors. Finally, different formulations of the energy equation are examined. 
A nonconservative form is found to be more accurate. 


I. Introduction 

I N the past tew years there has been a resurgence of interest in 
performing large-eddy simulations (LES)of flows of engineering 
interest. There are two roles for LES to play in the computation 
of such flows. First. LES can be used to test lower order models: 
k-t, algebraic stress, and full Reynolds stress models. LES can 
provide detailed data, which is difficult or impossible to measure 
experimentally and which is at much higher Reynolds numbers than 
can be reached by direct numerical simulation (DNS). Statistical 
data and physical insight gained from these simulations can be used 
to evaluate and improve the lower order models. With this approach, 
however, the subgnd-scale (SGS) model used in the LES has to be 
validated in order to ensure that the LES data are correct. 

Second. LES can be used as an engineering tool rather than as a 
research tool. With the expected increases in computer capabilities 
in the near future, especially from the use of massively parallel 
computers, it may be feasible to perform LES of flows of engineering 
interest. LES will remain an expensive tool, but it will likely be the 
only means of accurately computing complex flows for which lower 
order turbulence models fail. 

Recently, there has been much interest in using the dynamic 
SGS model to perform LES. The dynamic model was first intro- 
duced by Germano et al. 1 and was extented for use in compressible 
flows by Moin et al. : Since then, further refinements to the model 
have been proposed. 3-6 The main advantage of the dynamic model 
over other SGS models used in the past is that it requires little 
prior experience with the type of flow being considered. The model 
(dynamically) adjusts to the flow conditions by employing the re- 
solved large-scale information to predict the effects of the small 
scales. 

So far, the dynamic model has been mostly tested in incompress- 
ible turbulent flows and has been found to perform well. Moin et al. 2 
applied the dynamic model to compressible decaying isotropic tur- 
bulence and found that it performed well and better than using fixed 
values for the model constants. El Hady et al. 7 applied the model 


Received Jan. 5. 1995; presented as Paper 95-0355 at the AIAA 33rd 
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reserved. 
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to a transitional supersonic axisymmetnc boundary layer with satis- 
factory results. However, a number of issues regarding the use of the 
model in LES of compressible turbulence remain to be addressed, 
such as the ability of the model to capture compressibility effects 
without the need for explicit compressibility corrections. 

In addition, there are issues that need to be addressed in applying 
the dynamic model to inhomogeneous flows. The dynamic model 
requires filtering the resolved large-scale field. So far. it has been 
mostly implemented in turbulent flows that are homogeneous in at 
least two directions where the filtering can be performed efficiently 
(and exactly) in wave space using fast Fourier transforms. In more 
complex inhomogeneous flows, this is not possible, and some kind 
of discrete filtering has to be applied in physical space. In simu- 
lations of such flows, three-point explicit filters have been mostly 
used. A number of other spatial filters are available and require 
testing. 

The main objective of this paper is to examine the performance 
of the dynamic SGS model in the context of compressible decay- 
ing isotropic turbulence. The model is evaluated by making com- 
parisons with results from direct numerical simulations, as well as 
with reported “high 1 * Reynolds number, nearly incompressible ex- 
perimental data. The simulations are used to assess the capture of 
compressibility effects and to investigate issues regarding the im- 
plementation of the dynamic model for inhomogeneous flows. The 
reason for considering homogeneous turbulence is that the perfor- 
mance of the dynamic model can be evaluated separately from the 
effects of inhomogeneity. 

IL Mathematical Formulation 

A. Governing Equations 

In LES one computes the motion of the large-scale structures, 
while modeling the nonlinear interactions with the small scales. 
The governing equations for the large eddies in compressible flows 
are obtained after filtering the continuity, momentum, and energy 
equations and recasting in terms of Favre averages. The filtering 
operation (denoted by an overbar) maintains only the large scales 
and can be written in terms of a convolution integral. 

f(x)= I G(x—x')f(x‘) dr' (l) 

Jo 

where / is a turbulent field, G is some spatial filter (usually a sharp 
cutoff defined in Fourier space) of width equal to the grid spacing, 
and D is the flow domain. 



SPYROPOLTLOS AND B1-AISD6LL 


991 


The resulting equations of motion are as follows: 
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represents the resolved-scale stress tensor. The effects of the small 
scales are present in these equations through the SGS stress tensor 
and the SGS heat flux. 


T if = piu,u, - u,u ; ) 
q, = p{u, T - u, T) 


( 6 ) 

(7) 


respectively, and require modeling. A tilde is used to denote Favre 
averages if — p f j p). Also p is the density. T is the temperature, u, 
is the velocity vector, and k is the thermal conductivity. The specific 
heats at constant volume C v and at constant pressure C p are assumed 
m this study to be constant. The large-scale molecular viscosity jl 
is assumed to obey the power law whereas the 

large-scale pressure p is obtained from the filtered equation of state 
p = pRT . The molecular Prandtl number Pr is assumed to be 0.7. 
Note, that in deriving Eqs. (2—1), the viscous, pressure-dilatation 
and conduction terms were approximated in a similar fashion as by 
Erlebacher et al.* 
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The model requires the above averaging procedure, nn ad hoc 
solution, 1 to prevent numerical instabilities due to a simplification 
made in the derivation of the expressions for the model coefficients. 
For flows having a direction of homogeneity, spatial averaging is 
usually performed along that direction. Forthe case of homogeneous 
turbulence, this results in volume averaging and is the approach 
taken in this study. For this type of flow, the model coefficients vary 
onlv with time. 


B. Dynamic Modeling of the Subgrid Scales 

The SGS terms [Eqs. (6) and (7)) are modeled here using a com- 
pressible flow version of the dynamic SGS model of Germanoet al.‘: 
for details, see Refs. 2 and 3. The model involves three coefficients, 

C. C , . and Pr t . They are automatically adjusted, as time progresses, 
based on the resolved flowfield information with the aid of a second 
filter (test filter G) that has a filter width coarser than the gnd used 
to perform the computations. 

A refinement to the Moin et al. : model has been proposed by 
Lilly. 1 The two versions differ only on the type of a contraction used 
to determine uniquely the model coefficients, as is described later 
in this section. Both versions were tested in simulations, and com- 
parative results are presented in Sec. III.B. The results presented in 
other sections were obtained using the Lilly contraction. 

The model parameterization for the SGS stress and the SGS heat 
flux is given by 

r„ - ir u fty = -2CM J |5|(5„ - $$»**,) (8) 

. r»* =2C,pA : |5| 2 (9) 


<7. = 


pCA : |5l <77 

Pr, i)x, 


( 10 ) 


where 


C. Computer Implementation 

The numerical method for the direct and large eddy simulations 
employed a pseudo- spectral Fourier collocation scheme for spatial 
di fferencing and a third -order Runge-Kutta method for advancing i n 
time. 9 The validity of the numerical implementation of the dynamic 
model was established by performing a priori tests similar to those 
by Moin et al. : and comparing with their reported data. 

III. Results 

A. Capturing of Compressibility Effects 

The ability of the dynamic model to capture compressibility ef- 
fects was examined by performing LES of decaying isotropic turbu- 
lence and comparing with the results obtained from DNS. A number 
of simulations were conducted at different initial levels of com- 
pressibility and Reynolds numbers. The cases listed in Table l were 
considered. 

The level of compressibility of the initial fields was controlled 
by either varying the initial turbulent Mach number M, (cases 1-3). 
or the fraction of the turbulent kinetic energy initially contained 
in the dilatational velocity field x (cases 4—6). The effect of the 
turbulent Reynolds number Rer on compressibility was examined 
in cases 7. 8. and 3 and also in cases 4, 2. and 9. [Here. \f, = q/c . 
X = i q u !q) 1 . and Rer = pq l /(€p), where q is the rms magnitude 
of the fluctuation velocity, c is the mean speed of sound. q A is the 
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Table 1 Case parameters 


Case 

-Vf f 

X 

Rej 

1 

0.2 

0 

2742 

-> 

0.4 

0 

2742 

3 

0.6 

0 

2742 

a 

0.4 

0 

2157 

5 

0.4 

0.1 

2157 

6 

0.4 

0.2 

2157 

7 

0.6 

0 

735 

8 

0.6 

0 

2156 

9 

0.4' 

0 

6170 
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Fig. I Time evolutions of rms density fluctuations for cases 1-3; effects 
of initial M ( . 



Fig. 2 Time evolutions of rms density fluctuations for cases 4—6; effects 
of initial x- 





Fig. 3 Time evolutions of rms density fluctuations showing effects of 
initial Rer for cases a) 7, b) 8. and c) 3. 



Fig. 4 Time evolutions of rms density fluctuations showing effects of 
initial Rer for cases a) 4, b) 2, and c) 9. 

rms magnitude of the dilatational fluctuation velocity, and € is the 
dissipation rate of turbulent kinetic energy per unit mass.] 

All cases with purely solenoidal initial velocity fields ( x = 0) had 
uniform initial density and temperature fields, whereas the density 
and temperature fluctuations in the others were obtained from the 
isentropic relations and the condition for acoustic equilibrium (see 
Sarkar et ai. lft ). The initial three-dimensional energy spectrum for 
each case was of the form 

E(k)<xk9xp[-2(k/k l ,) 1 ] (20) 

where the wave number of the peak of the spectrum k n was set at 4. 
The LES were computed on (32) 3 grids, whereas the DNS were 
computed on (128) 3 grids. 

Good, and percentwise consistent, agreement in all statistics con- 
sidered was found between the LES and the DNS for all cases. This 
is shown, for example, in Figs. where the evolutions of the 
rms density fluctuations between the LES and the (filtered) DNS 
data for the above sets of cases are compared. (The time axes in 
these figures have been scaled with the initial eddy turnover time 
r, defined as the ratio of the lateral Taylor microscale and the rms 
fluctuation velocity in a direction.) Similar findings were obtained 
by comparing other statistical quantities, as well as one- and three- 
dimensional spectra, indicating that the dynamic SGS model seems 
to be capturing compressibility effects well for isotropic turbulence. 

B. Comparison of Two Model Versions 

In the preceding results, Lilly’s version 3 of the dynamic model 
was employed. The Moin et ai. 2 version was also tested for all 
-ses It predicted higher values for coefficient C, smaller values 
>r ^efficient C/, and similar values with Lilly’s version for Pr ,, 
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Fig. 5 Time evolutions of model coefficients showing effects of different 
contractions for case 6: a) C, b) C/ 9 and c) Pr t . 



Fig. 6 # Three-dimensional solenoidal and dilatations! energy spectra 
for case 6 at//r = 4.37; effects of different contractions. 

as is shown, for case 6. in Fig. 5. Overall, the Moin et al. model 
also performed well but provided higher amounts of dissipation 
than Lilly’s, as can be seen in the three-dimensional solenoidal and 
dilatational energy spectra for the same case shown in Fig. 6. This 
is most evident at higher wave numbers. The results are taken at a 
time when the turbuFent kinetic energy had decayed to one-fourth 
of its initial value. 

It should be noted that the values of Pr f obtained from either 
version of the dynamic model were about 0.4-0.6 when the initial • 
temperature field had a three-dimensional spectrum similar to the 
velocity’s. In contrast, when the temperature was initially set to 
be uniform. Pr t values higher than unity were predicted by the 
model. This behavior is believed to be due to differences in the 
initial transients of the temperature fields and is another indication 
for the need for dynamic modeling. 


C. Effects of Vmrrtftg the Tfetf FUtcr WWtfc 

The only adjustable parameter in the dynamic model is the ratio 
a = A/ A of the widths of the test and the gnd filter (see Sec. (I B). 
Based on a prion and a postenon tests of incompressible transitional 
and turbulent channel flow. Germane et al. 1 suggested a value of 
- tor future use. They also suggested further investigations using 
different types of flows. This value has since been adopted by other 
researchers and was used in most of the simulations presented here. 

The sensitivity of the results on the choice of the filter width 
ratio was also examined here for two cases of highly compress- 
ible isotropic decaying turbulence (cases 6 and 9 from Table 1). 
Five values of or were considered: 1.6. 16/9. 2. 16/7, and 8/3. 
These correspond to Fourier cutoff wave numbers for the test tiher 
of 10. 9, 8. 7. and 6, respectively. Note that the use of smaller 
or greater or values is undesirable, since it results in test-tiltered 
quantities that arc either almost unaffected by the filtering or contain 
only very large-scale information, respectively, and usually leads to 
ill-predicted model coefficients. Results for such cases are not pre- 
sented here. 

Noticeable differences in the evolutions of the model coefficients 
were found when different values for at were used in the simula- 
tions. as is shown for case 6 in Fig. 7. (The model coefficients were 
also calculated a priori from DNS results, and similar differences 
were seen.) 

Surprisingly enough, the choice of a seems to have only a very 
small effect on the LES results. 1 For instance, the rms density fluc- 
tuations tor case 6 (shown in Fig. 2 for or = 2) varied by less than 
for the five values of a considered {the differences are graphi- 
cally indistinguishable). This di fference was even smaller tor case 9. 
Best agreement with the DNS was obtained when ihis parameter was 
set to 2. 





width ratio for case 6: a) C, b) C/ t and c) Pr t . 
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Fig. 3 Thirt-dimensional temperature spectra from simulations of the 
Comte- Beilot and Corrsin experiment; efTects of modifying the convec- 
tive terms. 


D. Modification of the Convective Terms 

ft should be noted that the convective terms in the filtered momen- 
tum and filtered energy equations. Eqs. (3) and (4), were modified 
to the skew-symmemc form. 


3 l 3 1.3/1 3u; 


( 21 ) 


where / refers to pu, and pf , respectively. 

This modification reduces the effect of aliasing errors, 11 u which 
seems to be a bigger problem in LES than in DNS because the flow- 
fields are less well resolved. Shown in Fig. 8 are three-dimensional 
temperature spectra from LES of the (high Reynolds number) 
Comte-Bellot and Corrsin experiment 13 (see Sec. III.F.2 for a de- 
scription of the simulation). The LES without modifying the con- 
vective terms had a pile-up at high wave numbers, which led to 
instabilities, whereas the simulation with the modified terms was 
well behaved. 


E. Alternative Formulation for the Energy Equation 

In LES. there is also a choice tn the formulation of the energy 
equation. The most popular approach is the solution of a noncon- 
servative formulation (internal energy equation), since it requires 
only modeling of the SGS heat flux (see Sec. il.A). In contrast, if a 
conservative form uotai energy equation) is used instead, then mod- 
eling of additional SGS terms is required. An alternative procedure, 
which does not require any additional modeling, is to solve for the 
pseudo total energy, defined as 

<p = pCof + pu,Ui/ 2 (22) 

The resulting equation is 



This i)rm is not in strong conservation law form; however, we refer 
to it here as the conservative form. 

In LES of decaying compressible turbulence (case 6 from 
Table i) both formulations performed well, but the nonconserva- 
tive form was more accurate. The conservative form was tound to 
have somewhat of a pileup in the high wave number part of the 
density and temperature spectra, as is shown in Fig. 9. This is be- 
lieved'to be due to aliasing errors originating from evaluating the 
temperature from the pseudo total energy. 

F. EfTects of Employing Different Spatial Test Filters 

As was discussed in Sec. II. B, the dynamic model requires (test) 
filtering the resolved large-scale field. In flows that are homogeneous 
in at least two directions, the filtering can be performed efficiently 
(and exactly) in wave space using fast Founer transforms. In more 
complex inhomogeneous flows, this is not possible, and some kind 
of discrete filtering has to be applied in physical space. In the next 
two sections, the formulation and testing of such filters is presented. 



Fig. 9 Three-dimensional temperature spectra for case 6 at </r = 
2.13; efTects of reformulating the energy equation. 


l . Formulation of Spatial Filters 

Many different formulations of spatial filters were examined to 
determine their behavior in Fourier space. This was done because we 
believe it would be desirable for the spatial filter to have properties 
similar to the sharp cutoff filter. 

For uniformly spaced grids, the spatial filters are formulated as 
follows: 






4" f ) - n 
2 



n = 0 



(24) 


where f, is the function value at node j and /, is the corresponding 
filtered function value. N, is the number of implicit coefficients, 
which arc given by (a, ), and N f is the number of explicit coefficients, 
which arc (a„). 

The Fourier transform of the filtered function is given by f = 
QT \ where T is the Fourier transform of the function and Q is 
the filter transfer function. For the filter formulation (24). the filter 
transfer function is 



where k is the wavenumber and N is the number of grid points. 

The first filters to be considered are approximations to the top-hat 
filter in physical space. The top-hat filter is defined by 

I I if x - A/2 < x f < x 4- A/2 

(26) 

0 otherwise 

where the filter width A is assumed here equal to 2 Ax (as was 
suggested in Sec. III.C). The convolution integral in Eq. (i) can 
be approximated in various ways. Using three collocation points 
and the trapezoid rule gives ao = 1/2 and a\ — 1/2. With three 
collocation points one can improve the accuracy of the integration 
by using Simpson's rule, which gives an = 2/3 and a x — 1/3. If the 
integration is done analytically, then a spectrally accurate formula 
is obtained. 

The filter transfer functions for these schemes are shown in 
Fig. 10a against that of the sharp cutoff filter. The top-hat filter 
behaves very differently from the sharp cutoff filter and is more like 
the Gaussian filter in that some of the low wave number modes are 
reduced in amplitude whereas some of the high wave number modes 
are retained. 

A second set of filters is obtained by approximating the sharp 
cutoff filter in a least squares sense, with certain constraints imposed. 
Three constraints are considered. Constraint l is that the filter be 
mean preservering. This is important because the filter is supposed to 
separate large and small scales and should not alter a uniform field. 
For the filter to be mean preserving,- it is necessary that Q(0) = 1. 
Constraint 2 is that the end point of the filter transfer function is fixed 
to be zero; i.e.. Q(k — N/2) = 0. This constraint is not required 
on any physical basis; however, we believe it is a desirable property 
of the filter transfer function. Constraint 3 is that the ratio of the 
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Fig. 10 Comparison of filter transfer functions for various spatial filters: a) top-hat filters, b) sharp cutoff least squares filters, c) Lele filters, and 
d) compromise filter. 


filter width to the grid spacing, A/ Ax, be 2, since this is the chosen 
standard value. 

To implement constraint 3, one has to define the filter width for 
a given filter. We have chosen to define the filter width in a manner 
analogous to how the integral length scale is defined for a two-point 
correlation. For a discrete filter this corresponds to 

A = V 4- Ax (27) 

,-i G| 

Since the filters are mean preserving 

G j = Gn — l ( 28 ) 

) = t 


filter and the Gaussian filter; however, for the sharp cutoff filter such 
a definition gives an infinite filter width, and for some of the other 
filters considered in this study, an imaginary filter width is obtained. 
Therefore, the definition given in Eq. (31) was adopted. 

The least squares problem is formulated as minimizing the inte- 
grated square difference between the filter transfer function $'(*') 
and the transfer function for the sharp cutoff filter £' co (£')- ^ ree 
parameters are the filter coefficients \a n ) and (a m ). From constraint 
1 we have 



Constraint 2 gives 


and 



(29) 


This definition is dependent on the number of grid points N , which 
is undesirable. A more useful formulation can be obtained by letting 
k' = 2k/ N . A k' = 2A k/N (where A k = 1) and <?'(*') = 5(*). 
where Q{k) is given by Eq. (25), so that 



g'(k') A t£ 


(30) 


Taking the limit N-+ oo and using the symmetry of Q' results in 



S'(A') d*' 


(31) 


This definition is consistent with the length scale obtained from the 
cutoff wave number when a sharp cutoff filter is used, which is the 
reason it was chosen. An alternate definition in terms of the second 
moment of the filter function was used by Leonard 14 for the top-hat 


v c - i 

g '( k ' = 1) = 0 = ]T <2, (33) 

a =0 

These two equations can be solved for a^ t ^\ and av f - 1 in terms 
of the other parameters, so that the numbers of degrees of freedom 
is reduced by two. Constraint 3 was enforced by using a penalty 
function. The function to be minimized is 

E = jf [£'(*') - dk' + A^±-2^ (34) 

The Nelder-Mead simplex search algorithm 15 was implemented to 
solve the minimization problem. The integrals needed to evaluate E 
were calculated using Romberg integration with the error tolerance 
set to 1CT 14 . The constant A was set to 10 10 so that constraint 3 
was met with an error on the order of 10~ 6 . For explicit filters, 
the integrals can be evaluated analytically, and the minimization 
problem can be solved exactly. This was used as a check of the 
numerical procedure. For this case. Eq. (31) reduces to 

A/ Ax = l/ao 135) 
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Table 2 Summary of spatial Alton tested \m sfan utotio c» 


Filter 



<*2 

J 3 



1 3-pt explicit, trapezoid 

05 

0.5 





2 3-pt explicit. Simpson 

0 6666667 

0 3333333 





3 5-pt explicit, least squares 

0.4726761 

0 5819719 

-005464480 




4 7-pt explicit, least squares 

05 

06744132 

0 

-0.1744132 



5 5-pt implicU-7-pt explicit, least squares 

05 

0 8105146 

0.4830535 

0.1725389 

0 

09661071 

6 5-pt implicit-7-pt explicit. Lele 

05 

0 75 

0.3 

0.05 

0 

0.6 

7 5-pt i mpl icit— 7- pc explicit, compromise 

05 

0.8056734 

0.4684092 

0.1627358 

0 

0.9368185 



Fig. 11 Time evolutions of rms density fluctuations for case 6; effects 
of spatial Alters: a) explicit and b) implicit 

The resulting filter coefficients are given in Table 2. The three- 
point explicit filter is the same as the trapezoid filter. No five-point 
explicit filter was found that met all of the constraints, and so con- 
straints 2 and 3 were dropped for the five-point explicit filter shown. 
The highest order explicit filter considered is a seven-point filter. A 
number of implicit filters were found. The one discussed here is a 
five-point implicit-seven-point explicit filter. 

The least squares filters are compared to the sharp cutoff filter 
in Fig. 10b. As the number of filter coefficients is increased, bet- 
ter agreement with the sharp cutoff filter is obtained, as expected. 
However, the filter transfer functions are oscillatory and display 
large amplitude overshoots, which is believed to be undesirable. 

A third set of filters is obtained following Leie 16 in which a dif- 
ferent set of constraints are placed on the filter transfer function. 
To compare these filters to the sharp cutoff filter with A/ A* = 2, 
the extra constraint Q(k = jV/ 4) = 1/2 was imposed. The filters 
obtained are shown in Fig. 10c. In general, the filters are smooth, 
monotonically decreasing from l to 0. However, the filters are not 
as sharp as the five-point implicit method obtained from the least 
squares approach. Note that the three-point explicit filter is the same 
as the trapezoid filter. 

The filters obtained from Leie's formulation are very smooth, 
whereas the filters found from the least squares approach are much 
sharper but exhibit oscillations. In an attempt to reach a compromise 
between the two filters, we constructed filters obtained by combining 
the filter coefficients in the following way: 

a n = i/a 1 * 4* (1 - or„ = i/ot* + (l - ^)a^ ic (36) 

Figure lOd shows the transfer function for the five-point implicit- 
seven-point explicit filter with ^=0.92. This compromise filter 



Fig. 12 Three-dimensional energy spectra for case 6 at t/r - 4.37; 
effects of spatial Alters: a) explicit and b) implicit. 

gives a transfer function that is sharper than that from the Lele 
formulation but does not have as large an overshoot as that from the 
least squares approach. 

2. Testing of Spatial Filters 

A number of spatial filters from the preceding section were tested 
in simulations. TTie Alters considered are shown in Table 2 (see also 
Fig. 10). Homogeneous Aow test cases were chosen, so that the 
effects of these filters can be compared to that of the sharp cutoff 
filter, defined in wave space, as well as against available DNS or 
experimental data. 

Results from simulations of a highly compressible decaying 
isotropic flow (case 6 from Table 1 ) are presented first. The histones 
of the rms density fluctuations are compared in Fig. 1 1. Figure 12 
shows comparisons of three-dimensional energy spectra taken at a 
time when the turbulent kinetic energy had decayed to one- fourth 
of its initial value. The results from a DNS are also included in 
these figures for further comparison. The evolutions of the model 
coefficient C are shown in Fig. 13. All spatial filters are found to be 
more dissipative at higher modes than the sharp cutoff. All implicit 
filters performed well and, among them, filter 7 was found to be 
the best candidate. Most of the explicit filters were very dissipa- 
tive and performed poorly, with the exception of filter 4 which gave 
results very similar to those from the sharp cutoff filter. However, 
this filter predicted small negative values for the coefficient C f and, 
consequently, from Eq. (9), negative values for the subgnd-scale tur- 
bulent kinetic energy r M . This may, in some cases, lead to numerical 
instabilities, although it did not here. 

The spatial filters were also tested for the case of a nearly incom- 
pressible turbulence. The isotropic gnd-generated (high Reynolds 
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Fig. 13 Time evolution of model coefficient C for cmse 6; effects of 
spatial filters: a) explicit and b) implicit. 


Fig. 15 Three-dimensional energy spectra from simulations of Comte- 
Beilot and Corrsln's experiment at If/® /M * 98; effects of spatial filters: 
a) explldt and b) implicit. 



Fig. 14 Time evolution of turbulent kinetic energy from simulations of 
Comte-Bellot and Corrsin’s experiment; effects of spatial Alters: a) ex- 
plicit and b) implicit. 

number) turbulence experiment of Comte-Bellot and Corrsin 13 was 
simulated as a temporal decay on a (32) 3 grid. The initial veloc- 
ity fields for the simulations were purely soienoidal and designed 
to have the same three-dimensional spectrum as that reported at 
r (Jn/M = 42 in the experiment (where U n is the mean flow velocity 
and M is the gnd size). The initial pressure fields were computed 
from the incompressible Poisson equation, while the density was 
assumed to be uniform. The initial Mach number was set to 0.3. 



Fig. 16 Time evolution of model coefficient C from simulations of 
Comte-Bellot and Corrsin’s experiment; effects of spatial Alters: a) ex- 
plicit and b) implicit. 

The sharp cutoff filter was employed first in the dynamic model. 
The predicted time development of the turbulent kinetic energy and 
the three-dimensional energy spectra compare well with the avail- 
able experimental data, as shown in Figs. 14 and 15, respectively. 

The performance of the spatial filters was then examined. 
Figures 14 and 15 also present results from these simulations. The 
evolution of the model coefficient C predicted from the various spa- 
tial filters are compared against theone from the sharp cutoff filter in 
Fie. 16. Filter 6 was found to be the most dissipative of the implicit 
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filters. However, all of these filters performed well and similarly to 
the sharp cutoff filter. Filter 4, the seven-point explicit filter, again 
performed very well, in contrast to the three-point and five-point 
explicit filters, which provided great amounts of dissipation* 

For the Comte-Bellot and Comin simulation, we neglected C { on 
the basis that for nearly- incompressible flows c u is small compared 
to the thermodynamic pressure. 117 18 This was done to avoid insta- 
bilities in the computations that occurred during the initial transient 
period for some of the spatial filters. However, negligible difference 
on the results was observed when the sharp cutoff filter was used 
and C/ was computed from the dynamic model. 

IV. Conclusions 

Several issues involving the use of the dynamic SGS model in per- 
forming LES of compressible turbulent flows have been examined 
by employing the model in simulations and comparing with results 
from DNS or experiments. Decaying isotropic turbulence was con- 
sidered in order to evaluate the performance of the model separately 
from the effects of in homogeneity. 

We conducted a parametric study where the levels of compress- 
ibility of the initial flow fields were varied* The model, with its 
ability to adjust itself to the flow conditions, was found to predict 
well, from a statistical viewpoint, the bulk of the flow. The dynamic 
model was able to capture compressibility effects well and does not 
require any explicit compressibility corrections. 

In performing LES of inhomogeneous compressible turbulent 
flows using the dynamic model, the filtering operation required by 
the model can not be done in wave space buL, rather, has to be ap- 
proximated in physical space. We have examined the behavior of 
several implicit and explicit spatial filters in wave space. We also 
conducted simulations of highly compressible and high Reynolds 
number nearly incompressible cases of decaying isotropic turbu- 
lence to study the performance of such filters against that of the 
sharp cutoff filter defined conveniently ia Fourier space. The five- 
point implicit-seven-point explicit filters examined performed ex- 
tremely well and should be considered as good candidates for future 
use. However these filters would be more expensive to employ than 
explicit ones. The former filters have transfer functions that are 
much sharper than the other filters near the cutoff mode, but exhibit 
oscillations that as suggested by the simulations are not of seri- 
ous concern. The seven-point explicit filter gave good results also; 
however, it predicted small negative values for model coefficient 
Cj , which is undesirable. The three-point and five-point explicit fil- 
ters were very dissipative,' especially in the high Reynolds number 
simulation. 

Two versions of the dynamic model were employed and tested in 
LE£ The version by t-ifly^ was found tg & mere accurate than the 
version by MoiaetaJ. 1 ' ... T _ 

A number of simulations were also ggr fod nedtg study the effect 
of the filter width ratio— die only adjustable parameter in the dy- 
namic model. The results suggest ar greatlnsensitmty of the model 
to this parameter. Best agreem ent with the results from DNS was 
obtained when this parameter ^vas set to 2. 

Modification of the- ^onvectivg terms in the filtered momen- 
tum and filtered energy equations improves the accuracy of the 
simulations, as well as the stability of the numerical method. In 
addition, the nonconservative formulation, of the energy equation 
was found to be somewhat better than the conservative form. 
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Abstract 

The effect on aliasing errors of the formulation of nonlinear terms, such as the convective terms in the 
Navier-Stokes equations of fluid dynamics, is examined. A Fourier analysis shows that the skew-symmetric 
form of the convective term results in a reduced amplitude of the aliasing errors relative to the conservative and 
nonconservative forms. The three formulations of the convective term are tested for Burgers’ equation and in 
larae-eddy simulations of decaying compressible isotropic turbulence. The results for Burgers’ equation show 
that, while in certain cases the nonconservative form has the lowest error, the skew-symmetric form is the most 
robust. For the turbulence simulations, the skew-symmetric form gives the most accurate results, consistent with 
the error analysis. 


L Introduction 

Numerical simulation of nonlinear conservation laws, such as the Navier-Stokes equations of fluid 
dynamics, can give rise to nonlinear instabilities. This was pointed out long ago by Philips [14]. The 
instability was attributed to aliasing errors, in which the product of low order modes creates high fre- 
quency modes that when discretized appear as low frequency modes. This difficulty was overcome by 
Arakawa [21 for two-dimensional, incompressible problems using the vorticity-streamfunction formu- 
lation by developing a finite differencing scheme which satisfies certain global conservation properties. 
The modifications ensure conservation of kinetic energy and squared vorticity in the absence of ex- 
ternal and viscous forces and time differencing errors. It was found that by meeting these constraints 
the instability was suppressed. 

This idea was extended to three dimensional simulations by Kwak [11], who used primitive variables 
and rewrote the convective term in the momentum equation. 
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in the modified form 


1 JL 

2 3x 7 


\ UUj 


1 du, 

- [x, , 

2 ; dxj 


( 1 . 2 ) 


where use has been made of the conservation of mass equation for incompressible flow. duj/dx } = 
0. The formulations in (1.1) and (1.2) are referred to as the nonconservative and skew-symmetric 
formulations respectively. A third formulation is the conservative form. 


Using the skew-symmetric formulation in (1.2), Kwak showed that the integral conservation properties 
are maintained by the finite differencing scheme. Mansour, Moin, Reynolds and Ferziger [12] showed 
that the identity equatins (1.1) and (1.2) does not hold numerically and that this is the reason the the 
global conservation properties are not satisfied. One can show that the reason (1.1) is not equivalent 
to (1.2) numerically is the occurrence of aliasing errors. The skew-symmetric form of the convective 
terms were extended to the simulation of compressible flows in [4,8,15]. 

In [12] the rotational form of the convective term was considered. It is given by 



and is widely used in incompressible simulations. The rotational form ensures additional conservation 
properties including that of helicity [10]. Horiuti [10] performed numerical experiments of incompress- 
ible turbulent channel flow using a finite differencing scheme with the skew-symmetnc form and the 
rotational form, and concluded that the skew-symmetric form was superior. Horiuti attributed the poor 
performance of the rotational form to large truncation errors in the vicinity of the wall. Zang [16] per- 
formed spectral calculations of incompressible turbulent channel flow using the skew-symmetric form 
and the rotational form of the convective terms. He was also able to perform dealiased simulations in 
which aliasing errors were removed, and he showed that the poor performance of the rotational form 
was due to aliasing errors. 

Zang states that no theoretical analysis is available to explain why aliasing errors are reduced for 
the skew-symmetric form. In the current paper we examine the aliasing errors created by the nonlinear 
convective term and show, through a Fourier analysis, that the magnitude of the aliasing errors is 
reduced for the skew-symmetnc form compared to conservative and nonconservative forms. The 
analysis method for the rotational form becomes too complicated to draw any conclusions regarding 
it; however, the low aliasing errors for the skew-symmetric form can be understood. We then test 
the three formulations on Burgers’ equation and in large-eddy simulations of compressible isotropic 
turbulence. 

There are other methods of reducing aliasing errors besides reformulating the convective terms. 
Orszag [13] pointed out that aliasing errors for a product can be eliminated by truncating the Fourier 
coefficients of the velocity and the product using the so-called “2/3 rule." For some problems, such 
as simulation of compressible turbulence, which are best time advanced in physical space rather 
than in Fourier space, this procedure requires additional Fourier transforms and can be prohibitively 
expensive. An alternative which does not use Fourier transforms is due to Anderson [1]. With this 
procedure the velocity and the product are filtered in physical space using a digital filter. Aliasing 
errors are eliminated by forming the product on a refined grid. This procedure can also be expense e 
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if high accuracy ts required. Canuto et al. [5] point out that aliasing errors cease to be an issue if 
adequate resolution is used. However, for problems such as large-eddy simulation of high Reynolds 
number turbulence, in which the flow field is not well resolved, aliasing errors are important and 
increased resolution is not a possible solution. Therefore, the reformulation of the convective terms 
discussed in the current paper offers an alternative to more expensive dealiasing procedures. 


2. Error analysis 

We are interested in examining the aliasing errors that arise in using pseudospectral Fourier methods 
to solve nonlinear partial differential equations. The specific problem we are concerned with is the 
solution of the compressible Navier-Stokes equations, and so we take the nonlinear convective terms as 
the form of interest. In order to analyze the aliasing errors, we consider a mode! problem of evaluating 
a derivative of the form 



The convective terms in the momentum and energy equations (and in a passive scalar or chemical 
species equation) can be put in this form. 

Consider the product fg. Using a spectral Fourier method, the functions / and g are represented by 
a finite set of Founer modes evaluated at N discrete collocation points as follows 

/,= y; o.2) 

tz= — jV/2t , I 

!)r'= 9me' kmXt : 0 . 3 ) 

m=-.V/2+l 

where x } = (j — 1)£/;V, (j = 1, . . . , ,V), are the collocation points, k n = 2 ttti/ L are the wavenumbers, 
L is the length of the computational domain, and f n and g m are the Founer coefficients. The continuous 
functions f(x) and g(x) are approximated by a finite set of Founer modes, 

.v /2 

fix) = ^2 fne lknZ > {2A) 

n=-:V/2+l 
.V/2 

gix) — ^ dm^ kr 

m— — N/2 J r\ 

The product is 

,V'2 

f(z)<jU) = Y1 

n = -.V/2-r! rn 
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-s i 


1 2.6) 
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When this product is discretized, some modes will have values of (k n 4 k m ) which lie outside the 
range of resolved wavenumbers, which is given by 



for *V even. On the discrete collocation points these modes will be identical to modes that lie within 
the range of resolved wavenumbers. Thus the unresolved, high wavenumber modes are ‘'aliased” to 
resolved, lower wavenumber modes. (See [5] for a more complete description.) A Fourier mode which 
has a wavenumber outside the resolved wavenumber range is aliased to the mode within the resolved 
wavenumber range which has a wavenumber that differs from (Ar n 4* fc m ) by an integer multiple of 
2 ~N/L. We will denote the wavenumber of the resolved mode to which the mode {k n A-k m ) is aliased 
by ( k n 4 k m )° . A graph of (k n 4 k m )° as a function of (k n 4 k m ) is shown in Fig. 1. Note that if 
(jfc n -f k m ) lies within the range of resolved wavenumbers then the mode is not aliased to a lower 
frequency mode and (fc n 4 k m )° = (k n 4 A: m ). 

The discrete derivative is found by transforming the discrete product values, multiplying the Fourier 
coefficients by i = >/— T times the wavenumber, and inverting the transform. This procedure gives 


,V/2 ,V/2 

= Z E ' (kn + k m yf n g m j lk ^ k ’* ) ° x >. (2.7) 

For modes that contnbute to the aliasing error the Fourier coefficients are multiplied by the wavenumber 
of the resolved mode rather than the actual wavenumber. 
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It will be shown that the skew-symmetric form of the nonlinear term gives a reduced amplitude of 
the above aliasing error. Consider rewriting the derivative of the product as 


(fo) ->-(1 -*)( T-9 + 


dx 


dx 


( 2 . 8 ) 


Upon substituting the Fourier representation for f(x) and g(x) and discretizing, this becomes 


A/- A/_ 

£ £ i(o(A!n + fcn)° + (l + 

n = — AV-+1 m=— .V/2--1 


(2.9) 


where the fact that fc° = fc for the resolved wavenumbers has been used. Since, for modes that create 
aliasing errors, (A; n 4* Jk m )° # (fc n 4- /c m ), the value of the derivative depends on the parameter a. 
When evaluated analytically, (2.8) is of course equivalent to (2.1) and independent of a. However, the 
product rule for differentiation used to obtain this form does not hold numerically. 

Let k* = [a(Jfcn + fc m )° + (l -a)(fc n + &m)] be an effective wavenumber that multiplies the 
Fourier coefficients when a derivative is taken. For Fourier modes with (k n 4- k m ) within the range 
of resolved wavenumbers, there is no aliasing error, whereas for modes that create aliasing errors, k* 
will affect the magnitude of the aliasing errors. Therefore, it may be possible to set the parameter a 
to reduce the aliasing error. 

With a = 1 we recover the original conservative form of the nonlinear term, while with a = 1/2 
we have the skew-symmetric form, and with ot = 0 we have the nonconservative torm. A plot 
of k* as a function of (k n 4- k m ) tor the three values of a is shown in Fig. 2. (Only the range of 
(k n 4- fc m ) that gives rise to aliasing errors is considered.) As seen in Fig. 2, the skew-symmetric 
form gives values of k* that are small for modes with (k n A k m ) close to the range of resolved 
wavenumbers. Modes in this range come from a product of modes with intermediate wavenumbers 
or the product of a mode with a low wavenumber and a mode with a high wavenumber. The value 
of k* for the skew-symmetric form is larger for higher values of (k n 4- & m ), and these modes come 
from products of two high wavenumber modes. For problems of physical interest, the solutions have 
Fourier coefficients that decay at high wavenumbers, and aliasing errors coming from the modes 
close to the resolved wavenumber range generally have higher amplitude than aliasing errors coming 
from higher wavenumber modes. Therefore, it is more important to suppress the aliasing errors coming 
from the region close to the resolved wavenumber range. Since the skew-symmetric form has a smaller 
amplitude for the aliasing errors coming from these modes than the other formulations, it is expected 
to perform better. 

To determine if the reduction in aliasing errors is realized in practice we consider two problems. 
The first is the one-dimensional viscous Burgers’ equation, and the second is a large-eddv simulation 
of compressible isotropic turbulence. 

Although here we are considering the convective term, other nonlinear terms, such as diffusion 
terms with nonuniform transport coefficients, can be written in the skew-symmetric form. We believe 
aliasing errors would be reduced for these terms as well. 
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Fiii. 2. Effective wavenumber in the derivative of a product. Conservative form. skew-symmetric form. - 

nonconservative form. 


3. Burgers’ equation model problem 

As a model problem we consider the one-dimensional viscous Burgers’ equation, 

Ut -r UU X = UU XX1 < 

where a subscript denotes partial differentiation. Benton and Platzman [3] give several exact solutions 
and we choose their case ( 5 . 2 ) in which the boundary conditions are periodic and the initial conditions 
are given by u(z, 0 ) = — Rsin(x). 

Burgers’ equation was solved numerically using a Fourier pseudospectral method with a 5th order 
accurate Runge-Kutta-Fehlberg time advancement scheme [9]. The nonlinear convective term was 
formulated using the following three forms 


7 (tr) x - \ uu x , <3.3) 

uu r , ' < 3 ^ 

which are the conservative, skew-symmetric and nonconservative forms respectively. The specific case 
we examine has R = —20 and was integrated to t = 0.075. The time step was ci t = 0.0001 and 
the time differencing errors were checked to ensure that they were negligible compared to the spatial 
differencing errors. The calculations were done on a Cray YMP tor which the word length is 6-1 bits 
in an effort to keep round off error low. 

The numerical solutions for .V = 32 are compared to the exact solution in Fig. 3. and the maximum 
and r.m.s. errors are given in Table I. At first the results are somewhat surprising because, while the 
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Fig. 3. Solution to Burgers’ equation at t = 0.075 (N = 32). Exact solution. — , conservative form. skew-symmetric 

form. — • — . nonconservative form. 

Table 1 

Maximum and r.m.s. errors in the numerical solution to Burgers’ equation 




iV 

= 32 

N 

= 33 

Formulation 

a 

max 

r.m.s. 

max 

r.m.s. 

Conservative 

1 

Ml 

0.127 

1.10 

0,056 

Skew-symmetnc 

1 /2 

0.59 

0.056 

0.33 

0.018 

Nonconservative 

0 

0.46 

0.026 

3.14 

0.202 


skew-symmetric form performs better than the conservative form, it does worse than the nonconser- 
vative form. However, we found that if the calculations are done using N = 33 grid points, the results 
(given in Table 1) show the skew-symmetric form performs the best and the nonconservative form 
does the worst, by a large margin. This behavior seems perplexing; however, a detailed investigation 
gives some insight into the differences in the formulations. 

The convective term was examined in detail to see the differences in the formulations. The exact 
solution u(x. t) was calculated at t = 0.075 at the N collocation points, and then the convective term 
was formed using the three formulations. The magnitudes of the Fourier coefficients are shown in 
Fiss. 4 and 5 for N — 32 and N = 33 respectively. Also shown are the Fourier coefficients for the 
exact convective term (the Fourier coefficients were obtained by transforming the exact convective 
term on N = 256 collocation points) and a dealiased convective term, which is described below. 
From Fia. 4 one finds that for .V = 32 the Fourier coefficients from the nonconservative formulation 
lie the closest to those from the exact solution, while for N — 33 the Fourier coefficients from the 
skew-svmmetric form are the closest. (Because of the symmetry of the problem, the phase of the 
Fourier coefficients is ~~iZ. and all of the numerical solutions have the same phase of the Fourier 
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Fig. 4. Magnitude of the Fourier coefficients of the c; -ective term formed using the exact solution at t = 0.075 on .V = 32 
gnd points. Conservative form. *, skew-symmetric form. -h nonconservative form, o, dealiascd. — . exact, . 


isor- — -r- — ! 



Fia. 5. Magnitude of the Founer coefficients of the convective term formed using the exact solution at t = 0.075 on A - 33 
grid points. Conservative form. *. skew-symmetric form. -K nonconservative form, o, dcaliased, — . exact. 

coefficients as the exact solution. Therefore, examination of the magnitude of the Founer coefficients 
of the convective term gives a proper comparison.) 

For .V = 32 the maanitude of the Founer coefficients for the conservative formulation are larger 
than those of the skew-symmetric form and those for the nonconservative formulation are smaller, 
while for _V = 33 these relative positions are switched. The switching of the relative magnitude of the 
Founer coefficients in the numencal solutions for .V = 32 and N = 33 is due to the phase relationships 
among the Founer coefficients of the solution. To see this, consider the factor i k*f n g m in (2.9). The 
phased the Founer coefficients of the exact solution are rr/2 for even positive or odd negative mode 
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numbers, and — -/2 for odd positive or even negative mode numbers. So f n = isgn (n)( — \ )”j/ n j. For 
the modes that create aliasing errors, the factor can be written 

i[t*(fc n -f- k m y - (l - a)(k n + k m )] [i sgn(n)(- [ ) n \f n \ ] [i sgn(m)( - 1 )'"|</mlj 

= \[a(k n - k m y + (\ — a) (k n 4- A- m )] |/^] j<7m|( — I ) n *' n+l . (3.5) 

(Since the sign of n and m will be the same for modes contributing to the aliasing error. 
sgn(n)sgn(m) = l.) If the product of modes n and m are aliased to mode /, then l = n -f m-r sgn(H.V. 
The Founer coefficients of the exact convective term have phases that are ±~; 2 with the sign opposite 
that of the coefficients for the solution (with the exception of the l = ±[ modes). Pulling out a factor 
to account for the phase we have 

[isgn(i)(-I) (,+ l, ]sgn(0[o(Jfc ri + * m ) 0 + (l -a)(A: n (3.6) 

Therefore, the aliasing error will be multiplied by (± l ) depending on whether the number of collocation 
points, ;\\ is even or odd. The error is either added or subtracted depending on the the parity of .V 
and the sign of AT, which is opposite for a = 0 and a = I. 

As an example, consider l = 15. With N — 32 the Fourier modes that contribute to the aliasing 
error are (n T m) = {(-15, -2), (-14, -3), . . . , (-2, - 15)}. The factor in equation (3.6), without the 
phase term, is 

[a(15) + (l -a)(-l7)]|/ n ||g 7n !(-l)- 32 = [32a - 17] |/ a ||y m |. (3.71 

For a = l, 1/2. 0 the factor in the brackets is 15, -l and -17, which corresponds very well to 
the relative magnitudes of the differences between the Fourier coefficients of the numerical solutions 
and the dealiased result for mode 15 in Fig. 4. With N = 33 the Founer modes that contribute to 
the aliasing error are (n.m) - {(- 16, -2), (- 15. -3), . . . , (-2, - 16)}. The factor in equation (3.6), 
without the phase term, is 

[a(15) + (l — a)( — 18)]j /n| 0 = [ 1 8 — jja] |/„ ] \[],n I • (a -3) 

For a = 1, 1/2. 0 the factor in the brackets is -15, 1.5 and 18, which again correspond very well 
to the results for mode 15 in Fig. 5. Therefore, we can understand the large shift in the error in the 
convective term observed when N is changed. 

To see the aliasing errors produced by the three formulations they are compared to a dealiased result 
in which the aliasing errors have been eliminated. Forming the product au T creates Founer modes 

with mode numbers — <V, V, which would require 2 N collocation points to represent. A dealiased 

convective term is formed by starting with u on N collocation points, transforming to wave space, 
padding the Founer coefficients with zeros so that at least 2N Founer coefficients are used (we used 
256 modes, although this is unnecessary), forming the Fourier transform of u x , transforming u and u. 
back to physical space on the larger number of collocation points, forming the product, transforming to 
wave space, and keeping only the N Fourier coefficients. This is the method that was used to compute 
the curve labeled "dealiased convective term" in Figs. 4 and 5. One sees that the skew-symmetric form 
is very close to the dealiased form, as one would expect from the analysis in Section 2. Therefore, 
the skew-symmetric form of the convective term gives a viable alternative to dealiasing. Note that 
dealiasina is very expensive for the compressible Navier-Stokes equations because of the extra Founer 
transforms needed. Although dealiasing is used for incompressible turbulence simulations, it is not 
practical for the compressible equations. 
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However, as is clear in Figs. 4 and 5, the dealiased convective term is not the same as the exact 
convective term. This is because the exact solution cannot be represented by a finite number of Fourier 
modes, and so the Fourier coefficients computed using N = 32 or .V = 33 collocation points are not 
the exact Fourier coefficients. The difference is due to truncation error. Just as with the aliasing 
errors, the value of the truncation error is determined by the phase relationships of the exact Fourier 
coefficients of the solution. Because of the symmetry of the solution, the parity of .V can make a 
lanze difference in the truncation error, as is clear in the difference between the dealiased solutions in 
Figs. 4 and 5. With this understanding one can see that the low value of the error in the convective 
term for the nonconservative formulation with N = 32 is due to aliasing errors partially canceling 
truncation errors, whereas with A* = 33 aliasing errors add with the truncation errors to produce a 
very large error. 

From the above analysis we see that the low error for the nonconservative formulation on A r = 32 
grid points is fortuitous and is due to the particular symmetry properties of the solution. Although the 
skew-symmetric form of the convective term does not give the least error in all cases, it is the most 
robust. For problems, such as the simulation of turbulence, where the solution has a more random 
distribution of phases of the Founer coefficients than for the solution of Burgers' equation, the skew- 
svmmetric form is expected to perform better than the other formulations. In the next section we 
examine the behavior of the various formulations applied to large-eddy simulation of turbulence. The 
above example also points out the danger of drawing conclusions about numerical methods based on 
simple problems, although such problems are useful in understanding how a numerical method works. 


4. Large-eddy simulation 


The problem we are interested in is the simulation of turbulent Hows, specifically compressible 
homogeneous turbulence. One approach is direct numerical simulation (DNS), in which the time 
dependent Navier-Stokes equations are solved without any kind of turbulence model. However, because 
of the wide range of length scales that occur in turbulence, DNS is limited by current computer 
capabilitv to low Reynolds numbers. A means of simulating flows at higher Reynolds numbers is to 
perform large-eddy simulations (LES), in which the motion of the large eddies are resolved, while the 
effects of eddies smaller than the gnd spacing are modeled. Aliasing errors occur in both DNS and 
LES; however, since large-eddy simulations are by definition not well resolved, aliasing errors are a 
more important issue with LES than with DNS. 

The large-eddy simulations used are discussed in detail in [15], so only a brief overview will be 
siven here. The simulations are of decaying isotropic compressible turbulence. The computer program 
uses a Founer pseudospectral collocation method with a compact storage third order Runge-Kutta 
time advancement scheme. The computations were performed on a Cray C90. 

The convective terms in the momentum and internal energy equation were written in the general 
form 


a 4 (/,3j)+(l " a) ( /, l 


dxj 


9j 


where in the momentum equation f, = f>u, and = a,, while in the energy equation j, — i>L ■■ L 
and ijj = a,. (Note that the internal energy equation is solved in LES ot compressible turbulence. 



APVUtt 687 v. 960401 Printed : 15/04/ 1996 ; 14:32 


File : APNUM687. t«x; VTELX/P p. 11 


G A Blaisdell et a /. / Applied Numerical Mathematics 20 (1996) 1-13 1 | 



Fig. 6. Decomposed velocity spectra from LES of compressible isotropic turbulence. Solenoidal spectra. E s (fc), dilatational 
spectra. £*(Jfe). Conservative form. . skew-symmetric form, — , nonconservative form. - - DNS data, o. 



k 


Fig. 7. Temperature spectra from LES of compressible isotropic turbulence. Conservative form. , skew-symmetric 

form, — . nonconservative form, - * — , DNS data. o. 


rather than the total energy equation, because of the difficulty in modeling additional terms that arise 
with the use of the total energy. See [7] and [15] for a discussion.) Setting a = 1, 1/2 or 0 gives the 
conservative, skew-symmetric, and nonconservative formulations respectively. 

The conditions used are those of case (6) in [15] for which the initial turbulent Mach number is 0.4, 
and the grid has 32 3 points. In order to see the differences in the effect of the formulations, we examined 
the three-dimensional power spectral density of the velocity held, decomposed into solenoidal and 
dilatational parts, and the temperature field. E s (k), E d {k) and Er(k) respectively. The velocity spectra 
are shown in Fis. 6 and the temperature spectra are shown in Fig. 7, at a time when the turbulence 
has decayed for two initial eddy turnover times. Also shown for comparison are the results from a 
DNS usins 12S 3 grid points. The results for the conservative formulation and the nonconservative 
formulation show an unphvsical pile-up in the spectra at high wavenumber k. This behavior is due to 
aliasing errors. The skew-svmmetnc formulation resuits, on the other hand, show very little pile-up 
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in the spectrum. As an additional test, simulations were also done on a higher Reynolds number case, 
which corresponds to the experiment of Compte-Beilot and Corrsin [6], The conservative form and the 
nonconservative form resulted in numerical instabilities that halted the simulations, while the skew- 
svmmetnc form resulted in a simulation that was completed and compared well with the experimental 
data. These results show that the skew-symmetric formulation of the convective terms reduces aliasing 
errors and produces a more accurate solution than the conservative and nonconservative formulations. 


5. Conclusion 

The effect of the formulation of the nonlinear convective term on aliasing errors has been exam- 
ined. A Fourier analysis using an effective wavenumber shows that the skew-symmetric form of the 
convective term results in a reduced amplitude of the aliasing errors compared to the conservative 
and nonconservative forms. The three formulations of the convective term were tested for Burgers’ 
equation. It was found that in some cases the nonconservative form gives the lowest total error rather 
than the skew-symmetric form. However, a detailed Fourier analysis of the convective term shows that 
the behavior of the formulations is dependent on the phase relationships of the Fourier coefficients of 
the solution, and that the skew-symmetric form gives the most robust results of the three. The three 
formulations were tested in the large-eddy simulation of decaying compressible isotropic turbulence, 
and it was found that the skew-symmetric form gives the most accurate results, consistent with the 
error analysis. 
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Introduction 

In the past few years, there has been a resurgence of interest in performing large-eddy 
simulations (LES) of flows of engineering interest. There are two roles for LES to play in 
the computation of complex turbulent flows. First. LES can be used to study the physics of 
turbulence at higher Reynolds numbers than can currently be achieved with direct numerical 
simulation (DNS) and LES can aid in the testing and improvement of lower order engineering 
turbulence models. Second, it is hoped that LES can be used in the near future as an 
engineering tool rather than as a research tool. Although it will remain an expensive tool, 
it might be the only means of accurately computing complex flows for which lower order 
models fail. 

The majority of LES reported in the literature involve incompressible fluid flows that are 
homogeneous in at least two spatial directions. While computation of such flows has greatly 
contributed to the development of LES, the computation of more complex flows is required. 
In this paper, the method is applied to a spatially-developing compressible boundary layer 
flow. Several issues related to the effect of the numerical scheme on the simulations are 
investigated. A high-order, upwind biased, implicit, finite difference scheme is employed in 
the simulations and subgrid-scale (SGS) modeling is performed using the dynamic model. 


Mathematical Formulation 

In large-eddy simulation (LES) one computes explicitly only the motion of the large-scale 
structures. The nonlinear interactions with the small-scales are not resolved by the numerical 
grid and are modeled. The governing equations for the large eddies in compressible flows 
are obtained after filtering the continuity, momentum, and energy equations and recasting 
in terms of Favre averages. The filtering operation (denoted by an overbar) maintains only 
the large-scales and can be written in terms of a convolution integral. 


f(x 1 ,x 2 l x 3 ) = / n <?.■(*.■ 

JD t=i 




1 



where f is <1 turbulent field, (S, - is some spatial filter that operates in the z~th direction and 
has a filter width. A, (usually equal to the computational grid spacing in that direction) and 
D is the flow domain. 

The resulting equations of motion for the large eddies are as follows: 


dp , d(pZi) _ n 

dt dx { 


at 

dp i d (pit,) 
dt 


d{piii) d(pu,uj) 


dxj 


dp d 


dxi dxj 


_ ( dui du 3 2 duk . A 

~ ni 


_dui . ( dui du 3 2 35* , A dv^ d_ ( rdf _ 
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(2) 

(3) 

(4) 


The effects of the small-scales are present in the above equations through the SGS stress 
tensor and the SGS heat flux, 

Tij = p {upU : - U,Uj) , ( 5 ) 

q t = p (uiT — Uif) , (6) 


respectively, and require modeling. A tilde is used to denote Favre aveiages (/ pj j p ) • 

Also, p is the density, T is the temperature, u t is the velocity component in the i-direction 
and k is the thermal conductivity. The specific heats at constant volume, C v , and at constant 
pressure, C p , are assumed in this study to be constant. The large-scale molecular viscosity, 
p, is assumed to obey Sutherland s law, 


p _ = (t A 3 / 2 r 0 + f e 

po~[fo) T + T e ’ 


with a Sutherland constant f e = 198.6 ° R. The large-scale pressure, p, is obtained from the 
filtered equation of state, 

p=pRT. (8) 


The molecular Prandtl number. Pr, is assumed to be 0.718. Note, that in deriving Eqs. (Si- 
ll), the viscous, pressure-dilatation and conduction terms were approximated in a similar 
fashion as by Erlebacher et aid For example, the pressure dilatation term is approximated 
as follows: 




(9) 


In the above, the small-scale temperature dilatation terms in the parentheses are neglected, 
since they are expected to have a small influence on the large-scales compared to the SGS 
heat flux. q>. and also because they are very difficult to model. 2 - 1 - 3 


Subgrid-Scale Modeling 

The dynamic SGS modeling concept was introduced by Germano et aid for LES of 
incompressible flows and has attracted a lot of attention in the LES community during the 



recent years. Moin et al? extented the dynamic model to compressible flows and Lilly 3 
suggested a refinement to both models that is now largely employed. Since then, further 
refinements to the model have been proposed (Wong, 6 Ghosal et al.,‘ Piomelli et al. s ). 

The model for the deviatoric and isotropic parts of the SGS stress tensor is based on 
Smagorinski’s 9 and Yoshizawa’s 10 eddy-viscosity models, respectively. The model constants, 
however, are allowed to vary in space and time, and are computed dynamically, as the 
simulation progresses, from the energy content of the smallest of the resolved large-scales. 
This approach of calculating the model constants has been found to substantially improve the 
accuracy and robustness of the LES method, since the model constants adjust dynamically 
to the local structure of the flow and do not have to be specified a priori. In addition, 
it has been found from incompressible flow simulations, that the dynamic model provides 
the correct limiting behavior near solid boundaries, and adjusts properly by itself in the 
transitional or laminar regimes. Although it can not predict properly backscatter. it allows 
for some reverse energy cascade. A similar approach is followed for the SGS heat flux. 

Dvnamic modeling is accomplished with the aid of a second filter (referred to as the test 
filter'd) that has a filter width A, in the i - th direction, that is coarser than the grid used 
to perform the computations (A; > A,). 

The model parameterization for the SGS stress and the SGS heat flux is given by 


th ~ | T kk Sij = -2f.it (5, - -SkkSij^j • 

(10) 


T kk = 2C/pA 2 |5| 2 , 

(11) 


, df 

« = - k> di, ■ 

(12) 

where p t — CA 2 p|5|, k t — p t /Pr t , 

Si, = 0.5 (dui/dx, + diij/dxi), |5| = ( 

2 bijbij) , and 

A = (AxAyAz) 1 ^ 3 . 

The model coefficients are computed from 


C = 

<(A; - l&kSij) Mij) 

( MpqMpq ) 

(13) 

Ci = 

(A-i) 

(14) 

{'’AWT - 2A'/|Y) 
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l/3 /v 

where ~ denotes test-filtered quantities, A = (AjA?-^) (A, is the width ot the test filter 

in the /th direction), < > denotes averaging over the homogeneous spanwise direction, and 


= P'UiUj - PUj ■ 

p 


3 



Mi 


ij 
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(17) 

( 18 ) 


Ki = JSt - l-pui pT , (!9) 

P 

In the simulations, negative values for the eddy viscosity, p u and eddy conductivity. k t . 
were allowed, as long as the total viscosity (/it = ft + f*t) an< l tbe tota l th erma l conductivity 
<k T = k + k t ) were non-negative. In terms, this restricts the amount of energy back-scatter 
allowed, but avoids numerical instabilities due to anti-dissipation. A three-point top-hat 
filter (derived using the trapezoidal integration rule) was employed for the test filtering. 


Numerical Method 

The DNS code of Rai et al. n was modified to perform LES using the dynamic model. The 
DNS code solves the Navier-Stokes equations in non-conservative form. Since the problem 
considered here does not exhibit discontinuities, the same approach was taken m the LES to 
solve Eqs. (2)-(4). Spatial derivatives are computed using fifth-order- accurate upwind-biased 
differences for the convective terms and fourth-order-accurate central differences for the 
viscous terms. Fourth-order-accurate central differences are also used to compute the spatial 
derivatives in the dynamic model. Time advancement is performed using an iterative fully 
implicit second-order-accurate scheme. 12 Such schemes are unconditionally stable and allow 
for accurate advancement using much larger time steps than explicit schemes. However, they 
are more CPU intensive, since they involve the solution of a system of algebraic equations. 
In addition, upwind schemes are much more stable than central difference schemes, since 
they provide implicitly some artificial dissipation (this controls also aliasing errors). 


Reference Ca.se and Calculation Set-up 

The experimental configuration of Shutts et al . 13 was chosen as a test case. It is that of 
a zero-pressure gradient, flat-plate boundary layer flow at M = 2.25. The Reynolds nuinb ^ r 
based on inlet conditions is 635000/in. The adiabatic wall temperature is 580 °R and the 
temperature at the freestream is 305 0 R. 

The size of the computational domain and the type of boundary conditions are chosen 
the same as in the DNS. 11 for consistency. The computational domain is divided along the 
streamwise direction in three regions. The first region is 2.5in. long and contains the regions 
of blowing and suction, as well as transition. The second region has uniform spacing m x. 
is 'Nn long, and contains the turbulent region. The third region is 6in. long and gradua ly 
becomes verv coarse to artificially damp the turbulent fluctuations and ensure that the outlet 
boundary will be non-reflective. The domain is 0.35in. wide in the span, and 3in. tall along 

the wall-normal direction. .. 

Periodic blowing and suction is imposed on the flat plate, at a distance of O.om. 10 m 
the inlet, to trip the incoming flow to turbulence. The wall-normal velocity at the plate m 


this region is as follows: 


= Au.^f{z)g{:)h(t) . 
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where 

f(x) = \sin6(\ — cos9)j 27 1/,2 , 

Imax 

cj(z) — * ^ ^ Zisiu\2'kI{z j z max ^P/)) 

^mai 

MO = £ r - sm(27rm(/?i + <p m )) , 

m=0 

and / moi = 10, m max = 5, x a = 4.5in., = 5in., A = 0.04 is the disturbance amplitude. 

= 75000Hz, are random numbers (between 0 and 1), and z maT = 0.35in. 

A no-slip boundary condition, together with an adiabatic wall temperature condition, is 
imposed on the rest of the flat plate. The conditions at the inflow and outflow boundaries are 
supersonic, except for the subsonic portion of the boundary layer. At the inflow boundary, 
the dependent variables are fixed based on results from a laminar boundary layer analysis. 

A non-reflecting boundary condition is specified at the outflow boundary, as was mentioned 
above. Periodic boundary conditions are used in the homogeneous spanwise direction. The 
computational domain is chosen such that it is long enough along the span to ensure that the 
flow is homogeneous in this direction. Finally, a symmetry boundary condition is imposed 
at the upper boundary, which is located well outside the boundary layer. 

Preliminary Results 

A number of simulations were conducted to examine several issues regarding the accuracy 
of the LES method in computing spatially evolving compressible boundary layer flows using 
finite difference schemes. The cases considered are summarized in Table 1. 

Shown in figure 1 is the variation of the computed skin- friction coefficient, C f . along the 
streamwise direction from simulations conducted at several different grid resolutions. The 
results are compared against the experimental data, the turbulent correlation of White and 
Christoff, 14 and the data from a marginally resolved DNS. 11 As the grid was refined the 
accuracy of the LES results improved. Case 2 uses 1/16 the number of grid points used in 
the DNS. This coarse grid simulation clearly under-predicts the skin friction, and appears to 
be at a much lower turbulence level close to the end of the well resolved region (x 8. Sin.) 
When the grid was refined (case 3) the results improved. The number of grid points along the 
spanwise direction was finally doubled (case 4) to capture better the large-scale structures. 
That simulation employed 1/3 the number of grid points used in the DNS. 

The variation of the Van Driest velocity (normalized by the shear stress at the wall) with 
the normal distance from the wall at x=8.8in. is shown in figure 2. Refining the grid again 
improved the agreement with the compressible law of the wall. Although the grid used in 
case 4 is only about a factor of 3 coarser than the DNS, there seems to still be room for 
improvement in accuracy. 

The skin-friction distribution from simulations performed at higher disturbance levels 
(cases 5) is shown in figure 3. As expected, the location of transition is moved further 
upstream as the disturbance amplitude is gradually increased. However, no difference on 
Cj are observed at x=8.Sin.. indicating that the flow there is fully turbulent. Therefore, 
the differences seen in the LES in figures 1 and 2. are not due to any end stage (by-pass) 
transition phenomena. 


9 = 2k(x — x a )/{xt, — ,r a ) . 

I max 

, Z, = 1 , Z, = 1.25Z /+t , 

i-o 

TTlmax 

£ Tm = 1 , T m = 1.25J m+1 - 

m=l 


0 



The effect of the numerical method was examined by employing a lower order accurate 
scheme in the simulations. The convective terms were computed using third-order upwind 
differences. Second-order central differences were used in computing the diffusion terms 
and the derivatives in the dynamic model. A significant drop in the computed skin-friction 
coefficient was found when the lower order scheme was employed, as is shown for cases 7 and 
S in figure 4. This figure also shows that the lower order scheme required about 2.65 times 
more grid points to match the results of the higher-order scheme. 

The final paper will contain comparisons of other boundary-layer statistics, such as pro- 
files of turbulence intensities. In brief, they have been found to compare similarly to the 
above results. 

Overall, the poor performance of the LES is believed to be mainly due to the truncation 
errors from the upwind scheme, rather than due to the dynamic SGS model. These errors 
artificially damp the turbulence of the smaller resolved scales. Subgrid scales contain less 
energy than the grid scales. As a result, even accurate modeling of subgrid scales will not 
overcome the errors due to the finite difference scheme. Furthermore, since the dynamic 
model predicts the eddy viscosity and eddy conductivity based on the turbulence level of the 
smallest resolved scales, it provides insufficient amounts of turbulent transport. Since the 
highly accurate spectral methods are not appropriate for use in complex flows, a possible 
solution to the problem would be to maintain only the information on the grid scales that 
are accurately resolved by a finite difference type scheme, while modeling the effects of the 
scales that are omited (including, of course, the subgrid scales). This approach, however, 
will substantially increase the cost of the simulation, since it would require explicit filtering 
of the contaminated modes at each time step using a high-order digital filter, and the use of 
finer grids to ensure that the remaining resolved scales adequately represent the large-eddies. 


Conclusions 

A number of issues involved in the LES of a spatially evolving compressible boundary 
lavers are examined by conducting simulations using a high-order-accurate finite difference 
scheme and the dynamic SGS model. The computational grid was refined successively to 
improve the agreement of the computed turbulence statistics with the available experimental 
data and results from a marginally resolved DNS. The grids used in the LES were from 16 
up to 3 times coarser than the grid used in the DNS. The computational domain was found 
to be long enough for the flow to reach a fullv-turbulent state and free of transients due 
to the periodic blowing and suction mechanism employed to by-pass the natural transition 
process. 

The results suggest that the finite difference scheme has a direct effect on the effectiveness 
of the SGS model influencing greatly the accuracy of the simulations. The use of higher- 
order scheme is found to improve substantially the results, since it improves the capture of 
the smaller resolved scales. Furthermore, it is recommended to apply the model not only 
at the subgrid-scales, but also at the scales that are not properly resolved by the numerical 
scheme. 
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Table 1: Case parameters 
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Figure 1: Streamwise evolution of the skin friction coefficient: effect of grid size 
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Figure 2: Profiles of Van Driest velocity normalized by wall-shear velocity at x=8.8in.; effect 
of grid size. 
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wise evolution of the skin friction coefficient; effect of disturbance amplitude. 
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